# Solve the partial differential equation using Crank-Nicolson method.

34 ビュー (過去 30 日間)
Hana Bachi 2022 年 2 月 18 日

Hello everyone,
you may see my questions here from time to time. I'm a beginner in Mathlab and trying to solve some equations, then ask for advice if I accounter any errors.
I tried to solve a PDE using Crank-Nicolson method. I wrote a code but I'm still getting some errors, which didn't allow me to get the final results. any advice, please.
I attached the code and the question.
Thank you

サインインしてコメントする。

### 採用された回答

Abolfazl Chaman Motlagh 2022 年 2 月 18 日
Hi.
you forget to put index in x vector in line 70, in equation you calculate Ur. that's the error.
Ur(i,j) = 0.5*erfc(0.5*(x(i)-t(j+1))./(B*sqrt(t(j+1))))+0.5*exp(x(i)/B^2).*erfc(0.5*(x(i)+t(j+1))/(B*sqrt(t(j+1))));
and also you put j+1 in t, but j is from 1 to M=1651. but at j=1651 the t(1652) doesn't exist. so correct that too.
##### 7 件のコメント表示非表示 6 件の古いコメント
Abolfazl Chaman Motlagh 2022 年 2 月 21 日
Hi hana.
i check your code, i couldn't figure it out exactly. there were a lot of ambiguities for me. specially variable d that you change it in every loop but never use it!
i tried to solve it myself. since my first try failed, the problem really got me :) . i literally solve the problem (discretization) more than 5 times from scratch. finally it works. the linear system is really sensitive to one little mistake in calculations.
i attached my code and also my formula for solution. hopes it is what you need.
here's final solution for 2 different discretization and the analytical approximation the pdf provide:

サインインしてコメントする。

### その他の回答 (1 件)

Hana Bachi 2022 年 2 月 22 日
Hello Abolfazi,
thank you so much for giving from your time to solve this problem, I deeply appreciate it.
The d term in my code is RHS in yours.
I found that I forget to multiply some terms by B^2, and I added dt/2 somewhere to, this is why it was showing big difference between the two solution.
But I think yours make more sence then the one I got
here is mine after writing the code

サインインしてコメントする。

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by