Solve the partial differential equation using Crank-Nicolson method.

39 views (last 30 days)
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

Accepted Answer

Abolfazl Chaman Motlagh
Abolfazl Chaman Motlagh on 18 Feb 2022
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 Comments
Abolfazl Chaman Motlagh
Abolfazl Chaman Motlagh on 21 Feb 2022
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:

Sign in to comment.

More Answers (1)

Hana Bachi
Hana Bachi on 22 Feb 2022
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!